%%%-------------------------------------------------------------------
%%% File    : p21.erl
%%% Author  : Plamen Dragozov <plamen at dragozov.com>
%%% Description : 
%%% Let d(n) be defined as the sum of proper divisors of n (numbers 
%%% less than n which divide evenly into n).
%%% If d(a) = b and d(b) = a, where a ≠ b, then a and b are an amicable 
%%% pair and each of a and b are called amicable numbers.
%%%
%%% For example, the proper divisors of 220 are 
%%% 1, 2, 4, 5, 10, 11, 20, 22, 44, 55 and 110; therefore d(220) = 284. 
%%% The proper divisors of 284 are 1, 2, 4, 71 and 142; so d(284) = 220.
%%%
%%% Evaluate the sum of all the amicable numbers under 10000.
%%%
%%% Created :  5 Dec 2008
%%%-------------------------------------------------------------------
-module(p21).

%% API
-compile(export_all).


%%====================================================================
%% API
%%====================================================================

%The prime factors of a number X:
% X = (p1**k1)*(p2**k2)..(pN**kN)
%Any divisor of X:
% Divisors = [(p1**a1)*(p2**a2)...*(pN**aN) for a1 in (0,..., k1), a2 in (0, ...,k2), ...,aN in (0, ...,kN)]
%The sum of all divisors of X:
% DivisorsSum = (p1**0 + p1**1+...+p1**k1)(p2**0 + p2**1 +p2**2+...+p2**k2)...(pN**0 + ...+pN**kN)

%%--------------------------------------------------------------------
%% Function: solution(N) -> int()
%% Description: The sum of all amicable numbers bellow N
%%--------------------------------------------------------------------
solution(N) ->
    scan(N - 1, ets:new(lookup, []), N, 0). 

%%====================================================================
%% Internal functions
%%====================================================================
% Iterate through all amicable numbers and sum them if they are bellow
% the maximum. Use a lookup table for saving calculations.
scan(0, _, _, Acc) ->
     Acc;
scan(N, LookupTable, Max, Acc) ->
    Prev = ets:lookup(LookupTable, N),
    case Prev of 
        [] -> 
            Sum = divisors_sum(N),
            DeSum = divisors_sum(Sum),
            case N =:= DeSum andalso N =/= Sum of
                true when Sum < Max -> 
                    ets:insert(LookupTable, [{N, Sum}, {Sum, N}]),
                    scan(N - 1, LookupTable, Max, Acc + N + Sum);
                true ->
                    ets:insert(LookupTable, [{N, Sum}, {Sum, N}]),
                    scan(N - 1, LookupTable, Max, Acc + N);
                _ -> scan(N - 1, LookupTable, Max, Acc)
            end;
        _ -> scan(N - 1, LookupTable, Max, Acc)
    end.

%Calculates the sum of all divisors of N (excluding N itself).
divisors_sum(N) ->
    {D, P} = divisor_powers(N),
    divisors_sum(D, P, 1) - N.

divisors_sum([], [], Acc) ->
    Acc;
divisors_sum([DivsH| DivsT], [PowersH| PowersT], Acc) ->
    divisors_sum(DivsT, PowersT, Acc*pow_sum(DivsH, PowersH)).

%Calculates the sum:
%X**0 + X**1 + X**2 + ...+X**N
pow_sum(X, N) ->
    pow_sum(X, 1, N, 1, 1).

pow_sum(X, Max, Max, Exp, Acc) ->
    Acc+ X*Exp;
pow_sum(X, Pow, Max, Exp, Acc) ->
    Next = Exp * X,
    pow_sum(X, Pow + 1, Max, Next, Acc + Next).

%Returns a list of all N's prime factors and their powers.
divisor_powers(N) ->
    {D, P} = divisor_powers(N, 2, 0, [], []),
    {lists:reverse(D), lists:reverse(P)}.

divisor_powers(N, I, Counter, Acc1, Acc2) when I > (N div 2)->
    case true of
        true when I =:= N -> {[N|Acc1], [Counter + 1 | Acc2]};
        true when Counter > 0 -> {[N|[I|Acc1]], [1|[Counter|Acc2]]};
        _ -> {[N|Acc1], [1|Acc2]}
    end;
divisor_powers(N, I, Counter, Acc1, Acc2) ->
    case true of
        true when (N rem I) =:= 0 ->
            divisor_powers(N div I, I, Counter + 1, Acc1, Acc2);
        true when Counter > 0 ->
            divisor_powers(N, I + 1, 0, [I|Acc1], [Counter|Acc2]);
        _ ->
            divisor_powers(N, I + 1, 0, Acc1, Acc2)
    end.
